clc,clear;
close all;
warning off;
%数据处理
% 导入数据
filename = 'newdata.xlsx';  
data1 = readtable(filename);
% 提取年份和省份数据
years = data1.x__;
pdata = data1(:, 3:end); 
nyears=length(years);
% 转换表格数据为矩阵
pd = table2array(pdata);
% 输入矩阵 X (实际就是pd)
X = pd(1:end-1, :);
% 输出矩阵 Y (计算年增长率)
Y = diff(pd) ./ pd(1:end-1, :);
%环境矩阵
Z=ones(33,3);
% 数据标准化：将 X 和 Y 矩阵中的数据归一化到 0-1 之间(对每组数据都同样操作)
P = (X - min(X(:))) / (max(X(:)) - min(X(:)));
Q = (Y - min(Y(:))) / (max(Y(:)) - min(Y(:)));
% Z1 = (Z - min(Z(:))) / (max(Z(:)) - min(Z(:)));
%计算3阶段DEA
% for k=1:nyears
%         P=[P1(k,:),P2(k,:)];
%         Q=[Q1(k,:),Q2(k,:)];
%         S=[Z1(k,:),Z2(k,:)];
%         % 计算 DEA 效率(各省)
%         n=size(P',1);m=size(P,1);s=size(Q,1);
%         A=[-P' Q'];
%         b=zeros(n,1);
%         LB=zeros(m+s,1);UB=[];
%         for i=1:n
%               f=[zeros(1,m)  -Q(:,i)'];
%               Aeq=[P(:,i)',zeros(1,s)];beq=1;
%               w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
%               E(i,i)=Q(:,i)'*w(m+1:m+s,i);
%         end
%         theta=diag(E)';%评价
%         
%         % SFA回归分析
%         sflts = cell(size(X, 2), 1);  % 用于保存每个输入变量的回归结果
%         for i = 1:size(X, 2)
%             slack = X(:, i) - theta(i) * X(:, i);  % 计算松弛变量
%             model = fitlm(S, slack);  % 进行回归分析
%             sflts{i} = model;  % 保存回归结果
%         end
%         % 提取 SFA 回归的结果 (系数)
%         beta = zeros(size(S, 2) + 1, size(X, 2));  
%         for i = 1:size(X, 2)
%             beta(:, i) = sflts{i}.Coefficients.Estimate;
%         end
%         % 调整输入数据 X
%         P0 = zeros(size(X));
%         for i = 1:size(X, 2)
%             adjust = S * beta(2:end, i) + beta(1, i);  % 计算调整值
%             P0(:, i) = X(:, i) - adjust;  % 调整后的输入数据
%         end
%         
%         % 使用调整后的输入数据进行DEA分析
%         n=size(P0',1);m=size(P0,1);s=size(Q,1);
%         A=[-P0' Q'];
%         b=zeros(n,1);
%         LB=zeros(m+s,1);UB=[];
%         for i=1:n
%               f=[zeros(1,m)  -Q(:,i)'];
%               Aeq=[P0(:,i)',zeros(1,s)];beq=1;
%               w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
%               E(i,i)=Q(:,i)'*w(m+1:m+s,i);
%         end
%         theta1(k,:)=diag(E)';%评价
% end

% 计算 DEA 效率(各省)
n=size(P',1);m=size(P,1);s=size(Q,1);
A=[-P' Q'];
b=zeros(n,1);
LB=zeros(m+s,1);UB=[];
for i=1:n
      f=[zeros(1,m)  -Q(:,i)'];
      Aeq=[P(:,i)',zeros(1,s)];beq=1;
      w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
      E(i,i)=Q(:,i)'*w(m+1:m+s,i);
end
theta=diag(E)';%评价

% SFA回归分析
sflts = cell(size(X, 2), 1);  % 用于保存每个输入变量的回归结果
for i = 1:size(X, 2)
    slack = X(:, i) - theta(i) * X(:, i);  % 计算松弛变量
    model = fitlm(Z, slack);  % 进行回归分析
    sflts{i} = model;  % 保存回归结果
end
% 提取 SFA 回归的结果 (系数)
beta = zeros(size(Z, 2) + 1, size(X, 2));  
for i = 1:size(X, 2)
    beta(:, i) = sflts{i}.Coefficients.Estimate;
end
% 调整输入数据 X
P0 = zeros(size(X));
for i = 1:size(X, 2)
    adjust = Z * beta(2:end, i) + beta(1, i);  % 计算调整值
    P0(:, i) = X(:, i) - adjust;  % 调整后的输入数据
end

% 使用调整后的输入数据进行DEA分析
n=size(P0',1);m=size(P0,1);s=size(Q,1);
A=[-P0' Q'];
b=zeros(n,1);
LB=zeros(m+s,1);UB=[];
for i=1:n
      f=[zeros(1,m)  -Q(:,i)'];
      Aeq=[P0(:,i)',zeros(1,s)];beq=1;
      w(:,i)=linprog(f,A,b,Aeq,beq,LB,UB);
      E(i,i)=Q(:,i)'*w(m+1:m+s,i);
end
theta1=diag(E)';%评价


% 输出调整后的效率
% disp('调整后的 DEA 效率：');
% disp(theta1);
% 绘制效率的箱线图
% figure;
% boxplot(theta1);
% title('调整后的DEA效率值');
% xlabel('地区');
% ylabel('效率值');



